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Abstract 

This paper is concerned with the further development of a new numerical method, the 
space-time solution element (STS) method, for solving conservation laws. The present work 
focuses on the two-dimensional, steady, incompressible Navier-Stokes equations. Using first 
an integral approach, and then a differential approach, the discrete flux conservation equa- 
tions presented in a recent paper are rederived. Here (i) a simpler method for determining 
the flux expressions at cell interfaces is given; (ii) a systematic and rigorous derivation of 
the conditions used to simulate the differential form of the governing conservation law(s) 
is provided; (iii) necessary and sufficient conditions for a discrete approximation to satisfy 
a conservation law in E2 are derived; and (iv) an estimate of the local truncation error is 
given. 

A specific scheme is then constructed for the solution of the thin airfoil boundary layer 
problem. Numerical results are presented which demonstrate the ability of the scheme to 
accurately resolve the developing boundary layer and wake regions using grids which are 
much coarser than those employed by other numerical methods. It is shown that ten 
cells in the cross-stream direction are sufficient to accurately resolve the developing airfoil 
boundary layer. 



A New Flux-Conserving Numerical Scheme for the 
Steady, Incompressible Navier-Stokes Equations 

I. Introduction 

This paper is concerned with the further development of a new numerical method* 
for solving conservation laws. 1-4 The differences between the current method and the 
traditional finite-difference, finite- volume, finite-element, and spectral methods have been 
previously described. 1-4 The key differences may be summarized as follows: the current 
method (i) provides for a unified treatment of space and time; (ii) represents the local 
discrete solution through a Taylor series approximation that identically satisfies both the 
integral and differential forms of the governing conservation law(s); (iii) balances fluxes at 
cell interfaces as an integral part of the numerical formulation; and (iv) evaluates fluxes 
at cell boundaries using exact functional expressions (to the order of accuracy of the local 
expansions). 

Specific numerical schemes 1-4 based on (i) - (iv) above have demonstrated the follow- 
ing properties: First, flux conservation is satisfied both locally and globally. Second, high 
accuracy is achieved without coupling the solution across numerous cells or grid points. 
Third, cells communicate only with their immediate neighbors, in much the same way that 
discrete regions of a real fluid interact. Fourth, the discrete dependent variables and their 
derivatives are all treated in a unified and consistent manner. Finally, the schemes them- 
selves are conceptually simple and lend themselves to straightforward implementation. 

In this paper, we are concerned with the continued development of a new flux- 
conserving numerical scheme for the steady Navier-Stokes equations. 3 Previous results 
have demonstrated the ability of this scheme to accurately resolve internal boundary layer 
flows on coarse, uniformly-spaced grids. 

One purpose of this paper is to extend the internal scheme presented in [3] to external 
flow fields. The specific application that we consider is incompressible, laminar flow past a 
thin airfoil. In spite of the significant differences between external and internal flows, the 
scheme we propose here is a straightforward extension of the previously presented internal 
flow scheme. Through comparisons with the Blasius solution, we show that the thin airfoil 
boundary layer can be accurately resolved on sparse, coarsely-spaced grids. 

Another purpose of the present work is to further illumine the flux conservation scheme 
presented in [3]. For simplicity, we concentrate in this paper on the incompressible Navier- 
Stokes equations. Using first an integral approach, and then a differential approach, we 
rederive the discrete flux conservation equations. Here we (i) present a simpler method for 
determining the flux expressions at cell interfaces; (ii) provide a systematic and rigorous 
derivation of the conditions used to simulate the differential form of the governing conser- 
vation law(s); (iii) derive necessary and sufficient conditions for a discrete approximation 

*The Space-Time Solution Element (STS) Method, also called The Method of Space-Time 
Conservation Element and Solution Element 
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to satisfy a conservation law in £ 2 ; and (iv) provide an estimate of the local truncation 
error. 

In the next section, we derive the discrete flux conservation equations for the incom- 
pressible Navier-Stokes equations. We then construct a specific scheme for the solution 
of the thin airfoil boundary layer problem, and conclude with a discussion of numerical 
results. 


II. Numerical Formulation 

A. Conservation Laws for the Navier-Stokes Equations 

We consider the two-dimensional, steady, incompressible Navier-Stokes equations in 
dimensionless form. We assume that the viscocity p is constant, and denote the density 
by p and the Reynolds number by ReL, where Re L = pU ^ L . The parameters L and Uoo 
refer to some reference length and velocity, respectively. 

Let x and y denote the horizontal and vertical coordinates, respectively, of a two- 
dimensional Euclidean space £ 2 - Denoting the horizontal velocity component by u , the 
vertical velocity component by v, and the static pressure by p, the governing equations for 
the conservation of mass and momentum may be written in Cartesian coordinates as 5 



du dv 

dx dy 

(2.1) 

9 ( 2 
& (u 

Q 

+ P ~ r xx ) + q^( uv ~ T *y ) = 0 

(2.2) 

3 , 

Q 

~ r X y) + -^(V 2 + P ~ Tyy) = 0 

(2.3) 


where 



2 . 0 du dv 

Txx ~ ZRe^ dx dy 

(2.4) 


1 du dv 

Txy Re L ^dy dx 

(2.5) 


2 . 0 dv du 

Tyy ~ZRe L { dy dx } 

(2.6) 


(Although or may be eliminated from t xx and r yy using (2.1), we retain both terms 
here to be consistent with the compressible formulation presented in [3].) 
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By applying the divergence theorem to equations 2.1 - 2 . 3 , they may be written in 


integral form as 


0 

II 

T-8 

$ 

t-ci 

£ 

( 2 . 7 ) 

h-x m ' — 0 

Js(V) 

(2.8) 

0 

II 

T -3 

s 

t-C 

S' 

w 

( 2 . 9 ) 

where S(V ) is the boundary of an arbitrary region V in E 2 , and ds is equal to da ft, where 
n is the outward unit normal to S(V) and da is the length of a surface element of *S( V). 
The flux current density vectors, h M , h XM , and h YM , corresponding to the conservation of 
mass, x-momentum, and y-momentum, respectively, are given by 

r def , . 

h M — (u,v) 

(2.10) 

h XM d = (u 2 + p - T xx ,UV - T X y) 

(2.11) 

7* dtf f 9 1 \ 

h YM (u U T X y, V + P Tyy ). 

(2.12) 


Equations 2.7 - 2.9 thus express physical conservation laws for the conservation of 
mass and momentum in an axbritraxy region V of £2. 


B. Discrete Flux Conservation Equations — Integral Formulation 

Let E 2 be discretized by a mesh with nonoverlapping rectangular regions. We assume 
constant spacing Aar and Ay in the x and y directions, respectively. (See Figure 1.) Each 
of the rectangular regions in the mesh will be referred to as both a conservation element 
and a solution element. A conservation element is a discrete region in E 2 over which the 
discrete analogue of the integral conservation laws ( 2 . 7 ) - ( 2 . 9 ) will be imposed. A solution 
element is a discrete region in E 2 in which a local Taylor series expansion is employed to 
represent the physical solution. In general, they need not refer to the same discrete region 
(See [2] or [4]). A conservation element will be denoted by CE(i,j) and a solution element 
by SE(i,j). The boundary of a conservation element will be denoted by S(CE(i, j)), and 
the cell center by (x;,t/j). 
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We then assume that the u and v velocity components and the static pressure p 
can each be represented locally on a solution element by a two-dimensional Taylor series 
expansion about the cell center (x t ,yj) as follows: 


- u 0i0 + tXi >0 (x-Xj) + u 0 l (y-yj ) 


(2.13) 


+ u 2 , 0 (x-x ,) 2 + u lA {x - Xi){y -yj) + u 0 , 2 (y - yy) 2 


v(x,y;i,j) 


dtf 


Uo,o + v Xi0 {x - Xi) + V 0tl (y - yj ) 


+ i> 2 ,o (x-Xj ) 2 + v 1(1 (x - Xi)(y - yj) + u 0i2 (y - y 7 ) 2 


(2.14) 


def 

p(x,y,ij) = Po, 0 + PiA x ~ x i) + Po.i(y - yj) (2.15) 

+ Pa. o(x-X ,) 2 + p ltl (x -Xj)(y-y,-) + PoAy-Vj) 2 - 

For clarity, the i,j subscripts have been omitted from the coefficients in the Taylor series 
expansions. These coefficients are the unknowns to be solved for. 

The Taylor series coefficients are related to the derivatives of the discrete dependent 
variables at the cell center by 


= u 

(2.16) 

= dujdx 

(2.17) 

= du/dy 

(2.18) 

= ^d 2 ufdx 2 

(2.19) 

= d 2 ujdxdy 

(2.20) 

= \ d 2 u/dy 2 , 

(2.21) 


and similarly for v and p. 

The discrete analogue to equations 2.7 - 2.9 in E 2 is then given by 


/ 

JS(CE(i,j)) 


h 



= 0 


( 2 . 22 ) 
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ftxM * ds — 0 

(2.23) 

0 ftyjy * ds — 0 

JS(CE(i,j)) ~ 

(2.24) 

where 


l de f { \ 

h M ~ («>H) 

(2.25) 

T de f / 2 , \ 

hxM = (tf+ p - Irr, UV - T xy ) 

(2.26) 

T* de f f 2 . 

i}ym — (y Z Z xyi V, d* P Zyy ) 

(2.27) 

and 


Ixi = (2 du/dx - dv/dy ) 

(2.28) 

r XJ , = (du/dy + du/dx) 

(2.29) 

Ivs = 3 Re L ( 2d ~/ dy du/dx). 

(2.30) 


Equations 2.22 - 2.24 axe a coupled system of integral conservation laws in which the 
fluxes h M ■ ds , h XM ■ ds , and h YM ■ ds are conserved by way of the discrete variables u,v 
and p. Each equation takes the form 

<f h-ds = 0, (2.31) 

where the second-order expansion ft is a function of u, v and p. Since the form of the 
integrals in (2.22) - (2.24) are identical, the integrations can be carried out by way of 

equation 2.31, where 


h = f (h\h y ) 


(2.32) 


and 
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(2.33) 


h x (x, y; i,j) = h x 00 + h x 0 (x-Xi) + K^y-yj) 

+ h x 0 (x-Xi) 2 + h x lA (x - Xi)(y -yj) + K 2 {y-yj ) 2 

h y (x,y,i,j) = f hi o + *».(*- *,•) + hUy-yj) (2.34) 

+ h\ 0 (x-xif + h\ x {x - Xi){y -yj) + h y 2 {y-yj) 2 . 

It is understood that each of the coefficients in (2.33) and (2.34) are Junctions of the 
discrete variables tt 0 , 0 , u 0 ,o, Po,o, Ui.o, «i.o, Pi.o, etc.. For example, when h corresponds to 
h XM , the term h* 0 corresponds to the constant term of the expression u 2 + p - t xx , and 
similarly for the other terms. Once the results are obtained in terms of h, it is a simple 
matter to obtain the corresponding results for h M , h XM , and h YM by way of equations 2.25 
- 2.27. 

The boundary S(CE(i,j)) of each conservation element is a simple closed curve in E 2 - 
Consequently, the surface integration required in equation 2.31 can be converted into a line 
integration form. 1 With ds = dc n , where n is the outward unit normal to S{C E(i, j)) 
and do is the length of a surface element in E 2 , we have 

rfs d = dy i — dx j, (2.35) 

and 

h • ds = — h y dx + h x dy = g • dr (2.36) 

where 

g d ~ f ( -h y ,h x ) (2.37) 

and 

dr d = dx i + dy j. (2.38) 

The line integration is taken to be positive in the counterclockwise sense. If we denote 
the vertices of an arbritrary conservation element CE(i,j) by P, Q, R, and S as shown in 
Figure 2, we have 

m h ■ ds = ® g- dr 

J. S(CE(i,j)) ~ JPQRSij ~ 

d U [, T(PQ ) + J(QR ) + J(RS ) + J(SP)]. . (2.39) 
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[J( PQ)]i,j denotes the flux of h through the line segment PQij, and similarly for 
J(RS), and j(SP). We then have (omitting i,j subscripts) 


Similarly, 


J(PQ ) d = J - h y dx + h x dy 

- r 

Jxi+& 


/ 

Jxi- 


4 A 1 

i T o 


h y dx 


— h y dx + h x dy 


A y 

with y = y } + 



def 

ryj + ^i- 



A x 

AQR) 

- h x dy 

4-v " 

with 

x = Xi — 

2 

J(WS ) 

def 

yx, + ^ 

— I h y dx 

1 Ax ~ 

with 

1 

II 

Ay 

2 



J X% 2 





def 

fVi + ^ L 



A x 

ASP) 

/ dy 

4 - 

with 

X = Xj + 

2 


Carrying out the line integrations in equations 2.40 - 2.43, one easily obtains 


Ax 3 


J(PQ) = hlo + ^ 


Ay fry . Ay fry I fry 

, * 9 U 0,1 ' u o,o 


J(QR ) = - 


Ay 3 . _ . r Ax 2 


r Ax‘ , Ax . _ . . 

12 ^0,2 ~ Ay [— h 2 0 — h 10 + h 0 ' ( 


J(RS) = - 


Ax 3 fAy 2 

" Ax 


4 fcj, - + M,] 


Ay 3 


* / \ y* /\ 7! ^ 

J(SP) = + Ay [—A?.. + + K. ■ 


J(QR), 

(2.40) 

(2.41) 

(2.42) 

(2.43) 

(2.44) 

(2.45) 

(2.46) 

(2.47) 
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By virtue of equations 2.31 and 2.39 we require that 


J(PQ ) + J(QR ) + J(RS ) + J(SP ) = 0. (2.48) 

Thus, we obtain the flux conservation constraint 

K> + fc ?,. = 0 - ( 2 - 49 ) 

— * 

Imposing this condition, we obtain the following expressions for the normalized flux of h 
across the boundaries of CE(i,j ): 


J{PQ) 

Ax 


Ax 2 

TT 


h y 

U 2,0 


+ 




+ 


h y 

U 0,0 


J(QR) 

Ay 


Ay 2 

12 


K> 


+ 


Ax 2 


Ko 


Ax 


<0 


+ 


Ko 


J[RS) 

Ax 


Ax 2 

~12~ 


h y 

2,0 


+ 




+ 


h y 

U 0,0 


JjSP) 

Ay 


Ay 2 

12 


K, 


+ 


Ax 2 


Ko 


+ 


Ax 


Ko 


+ 


K,o- 


(2.50) 


(2.51) 


(2.52) 


(2.53) 


Equations 2.50 - 2.53 are a third-order-accurate representation of the flux through the 
boundaries of CE(i,j). 

The flux expressions above may now be expressed in terms of the discrete dependent 
variables of the Navier-Stokes equations by way of equations 2.25 — 2.27. Corresponding 
to (2.25), we have 

h = h M (2.54) 


so that 


h x = u 
h y = v 

and the mass flux conservation constraint corresponding to (2.49) is 

u li0 + u 0l i = 0. 


(2.55a) 

(2.55b) 


(2.56) 
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The normalized flux expressions for h M are 


J(P<3)„ _ Ax* Ay* 

— s 12 ” M + 4 


Ay 

~2~ u i,o d* v 0 0 


(2.57) 


J(QR) m Ay 2 , Ax 2 Ax 

= ■ «0,2 + U 2,0 — “T-^l.O + U 0,0 

Ay 12 4 2 


(2.58) 


J{RS)m _ Ax 2 Ay 


Ay 


Ax 


12 


^2,0 “I” a ^0,2 + ri U l,0 *>0, 


(2.59) 


J(SP) M Ay 2 , Ax 2 , Ax 

V - — = -r^-Uo. 2 + — — u 2 ' 0 + — «i.o + «o,o- 

Ay 12 4 2 


(2.60) 


Corresponding to (2.26), we have 


h = h 


XM 


(2.61) 


so that 


h x — u 2 + p — t 3 


(2.62a) 


h ^ ” xi i) r jy . 


(2.62b) 


The x-momentum flux conservation constraint corresponding to (2.49) is 

1 2 

«i,o«o,o + u 0il v o a + p 10 - + Ui,i) + g (4 u 2>0 — v ltl )] = 0, (2.63) 


and the normalized flux expressions for h XM are 


J(PQ) 


XM 


Ax 

Ax 2 Ay 2 . v 

— (u 2 ,0 ^0,0 + V 2i o Uq ,0 + U l,0 U l,o) + ^ ( U 0,2 V 0,0 + V 0 ,2 U 0,0 + U 0,l V 0,l) 

Ax 

2 L 


(2.64) 


+ ^[v 0 ,oU OA ~ U 0 . 0 U 1>0 - -^-(2u 02 + t>i,i)j - F K,,+g + “o,o «o,o 
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_ J(QR)xm 

Ay 

Ay 2 . 2 . Ax 2 o \ 

— (2 tXo,2 “I” U 0 ,l Po,2J ^ (2^2,0 ^0,0 "b ^1,0 "h P?, 0/ 


(2.65) 


"1" 2 ^ V 0,0 ^0.1 ^0,0 U l,0 (2«0,2 + U l.l)] "I - Po,0 (2 1,0 V 0,l) "4* U 0,0 


JjRS), 

Ax 


( 2 . 66 ) 


^£.2 /\ 7^2 

12 ^ 2 ’° ^°’ 0 ^ 2,0 ^ 0,0 + U ltQ V lt0 ) + “ (^ 0,2 ^ 0,0 “ 1 “ ^ 0,2 ^ 0,0 “ I ” ^ 0,1 ^ 0 , l ) 

2 ^0,0 ^0,1 UO'OUl'O (^^0,2 + V l.i)j JR,C ^^ 0 ’ 1 ^l f °) U 0, 0^0,0 


J(SP) 


X M 


Ay 

Ay_%n.. .. , -.2 , _ ^ , Aa Wo ..2 

12 


(2.67) 


‘(2 IZq, 2 'Uq.O ~j“ W 0|1 ”1“ jPq , 2 ) + * ( 2 ^ 2,0 ^ 0,0 ^ 1,0 -p2 t o) 


^ ^0,0 "^0,1 Uo.oUl.o ^ (2ti 0 ,2 + Ul,l)j “1“ Po.o 377c ^1,0 ^0,l) + 


u„ 


(In applying (2.50) - (2.53) to the x-momentum equation, we have replaced h* 0 with —ho , i 
because it gives a simpler expression.) 

Corresponding to (2.27), we have 


h = Hk 


( 2 . 68 ) 


with 


h x = uv — r 


Lxy 


h y = y 2 + p — T 


yy- 


(2.69a) 

(2.69b) 


The y-momentum flux conservation constraint is 

1 2 

V 1 , 0 U 0,0 + V 0 ,l U 0,0 + Po,l ~ 1(2 v 2,0 + U lf i) + j(4v 0i j-u i| i)] = 0, (2.70) 


10 


and the normalized flux expressions for h Y M are 


J(PQ)ym 

Ax 

— (2t; 2|0 V 0t0 + V l : 0 + P2,o) "H (2 Uo, 2 U 0,0 + V 0> 1 "f” Po,2) 


(2.71) 


Ay 

2 


1 1 2 

Uo.oVl.o + V 0 , C U,, 0 — ^-(2u Ii0 +U l t )J + p .,0 — - -^-(2 v 0 ,, — U l 0 ) + V. 0 


J(QR) 


YM 


Ay 


(2.72) 


Ay 2 
12 
Ax 
’ ~2 


Ax 2 


/ \ JL, s 

(U 0> 2 W 0,0 "I" W 0,2 U 0,0 + ^0,1 U 0,l) + ^ ( V 2,0 U 0,0 4“ W 2, 0^0,0 "I" V l,0 U l,o) 

|u 0 , 0 fl,0 + V 0i0 U lr0 - (2 V 2 ,o +«!,l)j — ^-( U 0,l + U l.o) + U 0,0«0,C 


J(RS) 


YM 


Ax 

Ax 2 Ay 2 

— (2u,,o U 0 .0 + <0 + Pz,o) + -^—(2 Wo, 2«0,0 + t>o,l + Po, 2 ) 


(2.73) 


+ [«o,o Vi,o + u o,o u i,o — (2t? 2 ,o + «i,i)j + Po.o — 2 ~Re~^‘ V °' 1 Ul '°^ V “ ,<! 


J(SP) 


YM 


Ay 


(2.74) 


Ay 2 

12 


Ax 2 


(u 0 , 2^0,0 + U 0 , 2 U 0 ,0 + Uo.iWo.i) + — (v 2 ,o«o,o + «2,oUo,o + U li 0 U li0 ) 


, 0 ^ 1,0 + t>o,o«i,o - ^-(2Vj,o+«i,i)] - jk“(“o,i +Vi.o) + u o, 0 ^ 0 , 0 - 

The following conditions are then satisfied on each conservation element: 


Ax ' 

+ — [«0 


J(PQ)m + J(QR)m + J(RS) m + J(SP) M = 0 (2.75) 

J(PQ)xm + J(QR)xm + J(RS)xm + J(SP)xm = 0 (2.76a) 

J(PQ) YM + J(QR)ym + J(RS)ym + J(SP)y M = 0. (2.76b) 
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Thus, the total mass and momentum flux out of each conservation element is exactly zero. 
Furthermore, the discrete variables u, v, and p satisfy the Navier-Stokes equations in int- 
egral form. 

The above formulation provides the framework through which local and global flux 
conservation is achieved. We may now turn our attention to the differential conservation 
laws (2.1) - (2.3). 


C. Discrete Flux Conservation Equations - Differential Formulation 


We begin with a consideration of the general conservation law 

VH = 0 . 


(2.77) 


Let 

h = (, h*,h ») (2.78) 

be defined and continuous on an open domain *D of E 2 ■ Suppose that the partial derivatives 
of h x and h y exist to all orders and are continuous on T>. Then, a necessary and sufficient 
condition for h to be a solution of (2.77) in V is that its partial derivatives satisfy 


d n h x 


+ 


d n h y 


= 0 


dx n ~^ dy^ dx n * 1 Sy*"*" 1 
for n = 1, 2, 3, ... and k = 0, 1,2, ...,n — 1. 


For n = 1 this gives 


and for n = 2 


dh x 

dx 


+ 


dh y 

dy 


= 0 


d 2 h x 

dx 2 dxdy 

cPh^ d 2 h y 

dxdy dy 2 


(2.79) 

(2.80) 

(2.81a) 

(2.81b) 


If h is a solution of (2.77), then (2.80) holds. Equation 2.79 follows by repeated differ- 
entiation of equation 2.80 with respect to x and y. The general result can be established 
by induction. 

Conversely, (2.80) shows that if h satisfies (2.79), then h is a solution of (2.77). 

On the basis of the above, we state the following theorem: 
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Theorem 2.1 Let h = (h x ,h y ) be defined and continuous on an open domain V of E 2 , 
and let the partial derivatives of h x and h y exist to all orders and be continuous on V. 
Then, h satisfies V • h = 0 throughout V if and only if its partial derivatives satisfy 


d n h x , , , d n h y 

dx n - k dy k ^ X,y ^ + ’ y 

for n = 1,2,3,... and k = 0, 1,2, ...,n — 1, 


for every (x,y) in V. 


Now suppose that h is a solution of (2.77), and that h x and h y are also analytic 
throughout V. Near any (x 0 ,y 0 ) in V, we have the uniformly convergent Taylor series 
expansions 


dh x dh x 

h x (x,y ) = h x (x 0 ,y 0 ) + —(x 0 ,y 0 )(x-x 0 ) + -^-(x 0 , y 0 ) (y ~ 2/o) + - 


dx 

d n h x 


(x-x 0 )" k (y — y 0 )* 


^ V^r O' "h* -mjc-XoT* 

- 2^ Z^[dx n - k dy k{Xo,yo) i (n-fc)! 


n=0 k = 0 


k\ 


(2.82) 


h y (x,y) = h y (x 0 ,y 0 ) + -^-(x 0 , y 0 ) (x - x 0 ) + -^-(x 0 ,y 0 ) (y - y 0 ) + ... 

r d n h y I (x - x 0 )^ (y ~ yp) fc 

= 2^ 2^[dx n - k dy k[Xo ' yo) \ (n-fc)! k\ 

n= 0 fe =0 * V 


(2.83) 


dh x dh x . . d 2 fc x w » . 5 2 /i z x . 

+ -^r(*-.»o)(*-*«) + ^(*.»»o)(y-y.) + - 


00 n — 1 


r d n h x si (x — X 0 ) n (y yo) 

= 2^ 2s[dx n - k dy k{Xo ’ yo) l (n-fc-1)! k\ 


n=l k = 0 


(2.84) 


dh y dh y . . d^h y . . . , \ , 

(x 0 ,y 0 ) + ^-«-(xo,yo)(s-*o) + 2-(®o,yo)(y-yo ) + ... 


aS 5 


dy “ dy y ~°'* QJ ' dxdy 

d n h y , j (x-x 0 ) n ~ fc (y-yo)^"^ 


d n hy , xUx-zor - (y - y 0 j v 
- Z_,2^[ax n - fc ay^ Io,yoJ J (n-fc)! (fc — 1)! 

^ ^ r 0"^ , /I (x - xo)”-*- 1 (y - y 0 )* 

- 2^ 2^[dx n - k -'dy k + l[X ° ,yo) \ (n-fc-1)! fc! 


(2.85) 


71=1 k — 0 
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, d° de f i a 1 d ±f a j a 1 d ±f a 
where dx ^ = 1, = 57, and ai o 3y i - gj- 


We now examine the effect of truncating (2.82) and (2.83). Let 

ix/ \ dnh * / 0 (* - x o) n_ * (y - y 0 )* 

" X,Z,[a x n-^ y fc( Xo ’^)J (n- k)\ k\ 

n=0 k=0 ^ V ' 


and 


( 2 . 86 ) 


l«/_ \ Y^T d n h y ( (x — x 0 ) n k (y — y 0 ) A 

hA x ’V) - (n-k)\ 


k dy l 


(n - Jfc)! A:! 


(2.87) 


n=0 k=0 

be the iVth-order Taylor series expansions of h x and /i y , respectively. Then 
A (i,y) h„(x,y) - (jy + i _ j.| 

where (zj, yj) is a point on the Une segment between (x 0 ,y 0 ) and (i,y). 6 Similarly, 

- w, t _ v^r v+'h’ 1 (x - ».)”+■-* (»-».)* 

^ (^,y) M*»y) - (N + i-k)i 


( 2 . 88 ) 


k = 0 


(N + l-k)\ k\ 


(2.89) 


Let 


d N+1 h x 


I dx N+l ~ k dy k 

Lood Nq of (x c 
where ( x,y ) is any point in JV 0 , then 


< M£ +ltk and 


d N+1 h y 


< K+k ,k 


dx N ^~*dy k 

in a neighborhood Nq of (x 0 , y 0 ) for k = 0, 1, ..., N +1. If Ax = (x— x 0 )andAy 

1 / \ • • . • XT 1 1 


JV + l 


Ax N+l-fc Ay fc 


= (y-y 0 ) 


(2.90) 


and 


N + l 


„ Ax JV+1-t Ay k 
I h y (x,y) — /i]J,(x,y)| < X 1 _ £)! k\ 


(2.91) 


Let Mat+i =sup{M* +1)fc ,M^ +l fc , fc = 0, 1, ...,iV+l}. Then we have the more conservative 
error estimates 


v s, . M n+1 (N + 2) [max(Ax, Ay)]^ 1 


(2.92) 
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and 


(2.93) 


,*(.,») - »«..*)! < M " l(JV + ^if x ’ Avr - ' 

If N is an even number, then [( ^^-)!] 2 is defined by 


K^)!] 2 = (t + 5) ! (t + !) ! 


=' & + + l + = (T 


(2.94) 


Equations 2.84 - 2.87 together with Theorem 2.1 show clearly that the function h N == 
(/i*,h5Jr) is a solution of the governing conservation law (2.77). Its order of accuracy is 
given by (2.92) - (2.93). It follows that a necessary condition for a local polynomial 
approximation h n to be the Taylor series approximation h n to h is that it be a solution 
of Equation 2.77. 

To proceed further, we need to formalize the notions of “a local polynomial approxi- 
mation” and “convergence.” 

Definition 2.1 Let (x 0 , y 0 ) be a point in E 2 , and let Ax and Ay be positive numbers. 
Let JV 0AlA y = j(x,y) : |x-x 0 | < ^ and |y-y 0 | < Then, a local discrete polynomial 

approximation is a function h N = (h x N ,h y N ) defined on NoAxAy by 


where h x n _ kk 


N n 


II 

E E h n-k,k( Al ’ A y) (* - *.)■"* (y - 

(2.95) 


ti= 0 k= 0 


h v N (x,y) = 

E E K-kA**, A « ) (* - *.)""* (y -».)*. 

(2.96) 


n=0 k=0 


and h y n _ k k are functions of Ax and Ay. 


Definition 2.2 A local discrete polynomial approximation h N converges to order K 
to the exact solution h of the conservation law V • Tt = 0 as Ax — ► 0, Ay — »• 0 if and 
only if for any e > 0, there exist numbers <5i > 0 and S 2 > 0 such that when Ax < Si and 


Ay < S 2 , 


^n-k,k(^ x ^y) 


d n h z 11 

dx*-*dy k ^ X °' yt ) (n — k)\ k\ 


(2.97) 


d n h y 1 1 

*?l < e (2 - 98) 

for n = 0,1,...,#, and k = 0,1,..., n, and the remainder R K = h N — h K -+ 0 as Ax — *■ 0, 


Ay -► 0. 
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We may now state and prove the two following important theorems. 

— # 

Theorem 2.2 Let e be any positive number, and let 0 < K < N. If hff(Ax, Ay) 
converges to order K to the exact solution h of the conservation law V ■ 'H = 0 as 

Ax —* 0, Ay — > 0, then for all sufficiently small Ax and Ay, 

|| h - h N (Ax, Ay)||oo < e. 


Proof: 


We have 


■4 -4 

* 

US* 

1 

* 

+ 

* 

1 

II 

(2.2a) 

^ || ^ ^JV || oo 

5: ||/i — ^n||oo 4“ || h N hjy || oo 

(2.2b) 

1 

* 

IA 

1 

* 

II OO + \\h K — h K || OQ + ||tf K — Rk 1 1 OO 

(2.2c) 


where R K = h N — h K and R K = h N — h K . The first two terms on the right hand side 

of (2.2c) are each less than | for sufficiently small Ax and Ay, and the third term is a 

— ♦ 

polynomial whose lowest order term is of degree K + 1. Since the coefficients of R K axe 
fixed, and R K —* 0, the third term is also less than | for sufficiently small Ax and Ay, 
and the theorem follows. 

Theorem 2.3 Let e be any positive number, and let 1 < K < N. If h^(Ax, Ay) 

— * — » — ¥ 

converges to order K to the exact solution h of the conservation law V • Tt = 0 as 

Ax — » 0, Ay — ► 0, then for all sufficiently small Ax and Ay, 

(n - *) h*_ M (Ax, Ay) + (k + 1) h y n _ k _ l k+1 (Ax, Ay) < e 

for n = 1,2, ... K , k = 0, 1, ...n — 1. 


Proof: 


Since h N ( Ax, Ay) converges to order K, for any n = 1, 2, ..., K, and k = 0, 1, ..., n — 1, 
we have ^ 

h n-k,k( Ax ^y) („_ it)! k\ 

d " hV ( 1 
• (X °’ y0) \(^k 




dx n ~ k 1 dy k+l 


(n — k — 1)! ( k + 1)! 


so that 


(n ~ k) h x n _ k>k { Ax, Ay) [g^kg^k (*•’ »•)] ( n _ k - 1)! A:! 


1 )! 


(2.3 a) 
(2.3b) 


(2.3c) 
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(fc + l)fcJ_ fc _i >Jk+ i(A*,Ay) -» m Q x n-k-iQyk+i ( x o’y<>)j ( n _ % _ i)i fc[ (^-^) 

=» (n-k)h* n _ k>k (Ax,Ay) + (k + 1) h s n _ k _ 1>k+l (Ax, Ay) {2.3c) 

r d n h x , , , d n h y , ,1 i 1 

[d^d^ {Xo ' Vo) + 0,y ° ; J (n — k — 1)! Ar!’ 

But since h is a solution of the conservation law, 



r d n h x d n h y 1 

(*•»/) 

by Theorem 2.1. Thus, 



(n-k)h’_ ltk (Ax,Ay) + (* + 1) h’ n _ k _, k+ ,(Ax,Ay) 0 

(Wj) 

as Ax — ► 

0, Ay — * 0, and the theorem follows. 


The meaning of Theorem 2.3 becomes clear when we consider the divergence of h N . 

With h N 

d =f ( h x N , h y N ) defined by (2.95) and (2.96), we have 



h x N {x,y) = ( n_ fc ) h n-k,ki x ~ x ^ n ~ k ~ 1( < y ~ y ^ k 

° X n= 1 fc=0 

(2.99a) 


■fi- h y N (x,y) = k ^n-fc,Jt( X_X °) n_fc (^ _ 2 /o ) fc_1 

° y n=l k=l 



= £ £ (* + 1) »*.*.,,*+,(* - x.)”- 1 - 1 (» - ».)* 

n=l fc=0 

(2.99b) 

so that 

d ^ 

V -h N = — h x N (x,y) + — h y N {x,y) = 

(2.100) 


EE [(» -fc-l.k+l ]C - (» - ».)*• 

n=l Jt=0 


According to Theorem 2.3, a necessary condition for h N to converge to order N is that 

(n + (t + 1) 1+ 1 - 0 (2.101a) 

as Ax, Ay — »■ 0, for n = 1,2, ..., iV, = 0, 1, n — 1. 

The implications of this are especially significant when it comes to numerical calcula- 
tions. In general, the mechanism whereby conditions (2.101a) are satisfied depends on the 
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particular numerical method being used. Finite-difference methods, for example, satisfy 
each of the conditions (2.101a) to a certain order through the difference approximations 
that are used. In this case, condition (2.101a) is satisfied to a given order, say order X, 
for n = 1. Then for higher values of n, conditions (2.101a) are satisfied to an order which 
is less than or equal to L. The higher order constraints expressed by (2.101a) do not 
result in independent conditions for a finite-difference scheme. Rather, the higher order 
constraints are automatically satisfied by virtue of the difference equations employed to 
satisfy (2.101a) corresponding to n = 1. 

On the other hand, when one solves for the unknown coefficients h^_ kk and h y n _ k k 
directly, as in the present approach, each constraint associated with (2.101a) represents 
an independent condition. Thus, to ensure that (2.101a) is always satisfied, one should 
require h x _ k k and h y n _ k k to satisfy 

(» - k) i;_ M + (t + 1 ) +1 = o (2.101b) 

for n = 1, 2, ..., N, k = 0, 1, ..., n - 1. That is, Kn should be a solution of the conservation 
law. As a result, conditions (2.101a) are not satisfied just to a certain order as in finite- 
difference methods, but rather are satisfied identically. 

When the coefficients h*_ k k and h y n _ k k are functions of intermediate variables (as 
in the case of the Navier-Stokes equations), each constraint associated with (2.101b) must 
be expressed in terms of the intermediate variables. We now show that it is possible to 
obtain these constraints directly without the need to re-express (2.101b) in terms of the 
intermediate variables. 

Let h = (h x , h y ) be defined on a domain V such that h x and h y are analytic at the 
point (x 0 , t/ 0 ) in X>, with series expansions that converge for all (x, y) in V. We then also 
have the convergent series expressed by (2.84) and (2.85). If h is a solution to (2.77), 
then its partial derivatives satisfy (2.79) for all (x,y), and in particular for (x 0 ,y 0 ). On 
the other hand, suppose that the partial derivatives of h x and h y satisfy (2.79) at (x 0 ,y 0 ). 
Since and are each convergent for all (x, y) in T>, so is their sum. Let RHS(2.84) 
and RHS(2.85) denote the right hand sides of (2.84) and (2.85), respectively. Then all the 
coefficients in the infinite series ^ = RHS(2.84) + RHS(2.85) are zero. Hence, 

2hl + ML = o, so that h is a solution of (2.77). We have established the following corollary 
to Theorem 2.1. 

Corollary 2.1 Let h = (h x ,h y ) be defined on a domain V of E 2 , and let h x and h y 
be analytic at the point (x 0 ,y 0 ) in V, with Taylor series expansions that converge for all 
(x, y) in V. Then, h satisfies V • h = 0 throughout V if and only if its partial derivatives 
satisfy condition (2.79) at the point (x 0 ,y 0 ). 

Since a polynomial is everywhere analytic, Corollary 2.1 applies to the expansions 
(2.95) and (2.96). To ensure that h N is a solution of (2.77), it is sufficient to require the 
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partial derivatives of h X \ and h y N to satisfy (2.79) at the point (r 0 ,y 0 ). 

We now illustrate this by way of the x-momentum equation. In view of Theorem 2.8, 
we require the second-order expansion h XM to be a solution of the governing conservation 
law. Applying Corollary 2.1 , we first require h XM to satisfy (2.80) at the cell center ( Xi,Vj ) 
of SE(i,j ): 

d , o . d , v 

■q~(w + P ~ In) + — I*y) — 


dx 


dy' 


2 4 (s) + - £<*■■> + [ l (s)ls + s[ l; (s)) - ^ = °- (2102) 


3y' 


We may then immediately write 

2 


2uo,o«i,o + Pi,o - + tt 0l i Uo.o + w 0>0 u 0j i - -^^-( 2 u 0 ,2 + v i.») — 


1 2 

Uj, o«o,o + Uo.i Uo.o + Pl.o — [(2u 0 , 2 + Vi.i) + 3 (4u 2 ,o — Vi,i)] = 0 (2.103) 


where the equality follows from (2.56). The first-order constraint expressed by (2.103) is 
identical to the i-momentum flux conservation constraint (2.63). 

We now require h XM to satisfy (2.81a) and (2.81b) at (xi,yj). Differentiating (2.102) 
with respect to x and setting the resulting constant term to zero, we obtain the second- 
order constraint 


2(2u 0 ,o u 2 , 0 + 0 + p 2 ,o) + 1 ^ 1,1 u o,o "b ^i,i u o,o -b u i,o v o,i -b u o,i t’l.o 0. (2.104) 


Similarly, differentiating (2.102) with respect to y, we obtain 

2(u 0 , oVo >2 + v 0 , 0 u 0(2 + Vo.iWo.x) + 2u u « 0 ,o + 2u, i0 u 0 ,i + Pi.i = 0. (2.105) 

In the same manner, we obtain the second-order constraints 

2n !iD -b v t>1 = 0 (2.106) 

2u 0 ,2 + u 1(1 = 0 (2.107) 

2(u 0> o^2,o + v 0 , 0 «2,o + v 10 Ui,o) + ’Zvy A v 0t0 + 2u li0 Vo.i + Pi.i = 0 (2.108) 


2 (2r 0 , 0 u 0 2 + Vg j + Po, 2 ) + Ui,i v o,o + v, tl u 0 0 + u 10 u 0 l + u o,i y t,o 0 (2.109) 

for h M and h YM , respectively. The first-order constraints for h M and h YM are identical to 
(2.56) and (2.70). 
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By requiring h M , h XM , and h YM to be solutions of the governing conservation law, 
we automatically ensure that local flux conservation is satisfied. This is a general result, 
which follows from the divergence theorem. Any function with zero divergence will satisfy 
local flux conservation. Thus, requiring the discrete flux vectors to obey the conservation 
law ensures the satisfaction of both local flux conservation and the necessary conditions 
for convergence to order TV. 

The formal order of accuracy of h M , h XM , and h YM as approximations to h M , h XM , 
and h YM is third-order. If third-order accuracy is maintained throughout the full develop- 
ment of a specific scheme, then the order of accuracy of the scheme remains third-order. 
An ideal error bound for the local expansions would then be given by (2.90) and (2.91). 


III. Application to the Thin Airfoil 
Boundary Laver Problem 

We now construct a specific scheme for the thin airfoil boundary layer problem. Con- 
sider the mesh shown in Figure 3. The airfoil lies on the x axis between x = 0 and x = 1, 
and the grid is stretched in the plus and minus y directions. (The nonuniform spacing 
does not introduce any new complication since the discrete equations presented in the 
previous section still apply if Ax and Ay are replaced by Ax,- and Ayy.) Note that our 
mesh includes both an upstream and wake region. 

Let N{ and Nj denote the number of solution elements in the x and y directions, 
respectively. There are six unknowns per cell for each of the three discrete variables u, v, 
and p. There are then 18 NjNj unknowns altogether, and we require lSTV^A^ conditions to 

have a closed system of equations. 

The first-order constraints (2.56), (2.63), and (2.70), together with the second-order 
constraints (2.104) - (2.109), immediately provide 9 NiNj conditions. These conditions 
ensure that the discrete flux vectors h M , h XM , and h YM are solutions of the governing 
conservation law (equation 2.77). 

We must also require that mass and momentum fluxes balance across each vertical and 
horizontal interface in the mesh. This gives 3 Nj(N{ — 1) + 3N,(N } — 1) — 3N a conditions, 
where N a is the number of solution elements between the airfoil leading and trailing edges. 

Boundary conditions account for an additional 47V, + ZNj + 47V a conditions. For each 
cell adjacent to the airfoil, we require the mass flux through the wall face, and the u velocity 
component at the midpoint of the wall face, to be zero. At the upstream boundary we 
specify the velocity, and at the downstream boundary we specify the pressure. Along the 
free-stream boundary cells, we specify zero y gradient conditions for u and v. 

Finally, we may set 

Pm = 0 
p 2 , o = 0 
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(3.1) 

(3.2) 


on physical grounds. These terms are negligibly small due to the nearly constant (zero) 
pressure gradient.* 

The above conditions ensure the satisfaction of local and global flux conservation, 
boundary conditions, and all other relevant physical requirements. We are now free to 
impose any other physically realistic condition to close the system. The number of condi- 
tions needed is N{(Nj — 1) — N a . This is precisely the number of horizontal interfaces in 
the mesh (minus those that coincide with the airfoil). Consequently, there is an additional 
degree of freedom in specifying horizontal interface conditions. Following [3], we require u 
to be continuous at the midpoint of each horizontal interface. 

By virtue of (3.1) and (3.2), there are 16 unknowns on each solution element. However, 
using the local constraints (2.56), (2.63), (2.70), (2.106), (2.107), and (2.109), six of the 
unknowns may be eliminated in terms of other variables. The total number of unknowns 
that must be solved for is then 10 N{Nj. 

The discrete boundary value problem outlined above is presented in the Appendix. 
The equations presented there are a coupled system of second-order polynomial equations 
in the unknown coefficients u„ „, u 0i0 > p 0) o> etc-- Solution of this system may be accom- 
plished very efficiently using Newton’s method. Because the thin airfoil velocity field has a 
preferred direction, the Newton iteration generally converges to the physical solution with- 
out difficulty. An initial guess of uniform flow is usually sufficient to ensure convergence. 
For flows with more complicated physics, it may be necessary to use a different solution 
technique (e.g., a time-iterative approach). 

The Jacobian matrix associated with equations A. 9 - A.29 can be arranged to have the 
structure shown in Figure 4. The form is the same as the Jacobian matrix associated with 
the internal flow scheme. 3 Because the matrix is nearly block diagonal, direct inversion is 
a suitable choice for the Newton iteration. By employing a sub-block pivoting strategy, 3 
the Jacobian matrix retains the same structure throughout the elimination process, and 
there is no matrix fill-in in the upper and lower diagonal blocks. 


IV. Numerical Results 

In this section we present numerical results from calculations of the thin airfoil flow 
field at a Reynolds number of 100,000. To assess the accuracy of our results, we compare 
with the Blasius solution. However, due to leading edge, trailing edge, and wake effects, 
the Blasius boundary layer profile will not compare well with a Navier-Stokes solution 
over the full length of the airfoil. On the other hand, at high Reynolds numbers, the 
Blasius solution does provide the correct boundary layer profile over a substantial part of 

*For many flows it will not be physically correct to impose conditions (3.1) and (3.2). The 
author is currently investigating ways of retaining both p M and p 2 „ for the more general 
case. 
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the airfoil. In particular, the Blasius solution may be used for comparison in the interior 
region where the leading and trailing edges are many boundary layer thicknesses away. 

Information from the Blasius solution may also be used in the construction of a mesh 
for numerical computations. If we let Re denote the Reynolds number based on the airfoil 
chord c, then the similarity relation 7,8 

” = -VE (4,1) 


becomes 



U/U c 


(4.2) 


/ Voo){x / c) 

where all qu an tities except rj and Re are dimensional. We may identify the free-stream 
velocity U with £/<», and take = 1 (since we assume constant viscocity). Equation 4.2 
then becomes 

Corresponding to any y, x and y are related through the equation 


(4.3) 


y 

c 



(4.4) 


Equation (4.4) allows us to estimate the location of the edge of the boundary layer for 
any fixed ^Ve take tj — 6.0, corresponding to ~ 0.999999, to be the edge of the 
boundary layer. An estimate of the boundary layer thickness S x at any — is then given by 



At the trailing edge we have the estimate 



(4.5) 


(4.6) 


For Re = 100,000, one gets 8 t . e . = 0.0268. 

In the construction of a mesh for numerical computations, one would expect that the 
free-stream mesh boundaries should extend to at least a distance of 8t.e. from the airfoil 
in order to obtain accurate results. In the present study, detailed numerical calculations 
were carried out with mesh boundaries located at ±y = 6 t . e ., l-56 t . e ., and 2.0 6 t . e ., where 
8 t e = 0.027. For each case, calculations were performed on a grid with uniform x spacing 
(Ax = .02) and nonuniform y spacing with exponential stretching. The y spacing at the 
wall (Ayu,) was determined on the basis of the boundary layer thickness near the leading 
edge of the airfoil. With Ax = .02, we should require the mesh spacing to be fine enough 
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to resolve the boundary layer at x — .01 (since we require the u velocity to be continuous 
at the midpoint of horizontal cell interfaces). The boundary layer thickness at x = .01 is 
approximately 0.00268 (based on equation 4.5). Because the present discretization repre- 
sents the u velocity by a quadratic polynomial, one would expect to be able to resolve the 
nearly linear boundary layer profile in the near wall region with only one or two cells. Using 
two cells, one obtains the estimate A y w = .00134. However, convergence problems were 
encountered by the Newton’s method for this value of A y w (with uniform flow conditions 
as the starting solution). We used the slightly larger value of Ay w = .0015 instead. 

Once the wall spacing and mesh boundaries are fixed, grid stretching can be used to 
reduce the number of cells required to reach the outer boundary. Since we compute the 
flow field on both sides of the airfoil, the number of cells available to resolve the boundary 
layer is half the total number of cells in the y direction. A systematic series of calculations 
performed with mesh boundaries located at ±y = <5t. e ., 1-5 a-nd 2.0 St.e. showed that 
no accuracy was gained by extending the mesh boundaries farther than iy = S t e from 
the airfoil. Consequently, the present discussion will deal only with results obtained from 
the smaller mesh. 

The first series of calculations were performed on grids with uniform x spacing of Ax = 
.02. The upstream boundary was located at x = —.12, and the downstream boundary at 
x = 1.5. The spacing was varied exponentially in the y direction, using as few as 16 cells 
with exponential stretching of 22.4%, and as many as 28 cells with exponential stretching 
of 3.8%. In general, the agreement with the Blasius solution improved as the mesh was 
refined. In Figure 5 we show two of the grids that were used. 

Figures 6-8 present comparisons of the discrete velocity u with the Blasius solution 
at various x locations along the airfoil. The numerical results in these figures were obtained 
from the grids shown in Figure 5. Note that at each x,- (the x nodal points), u(i,j ) = 
[it 0>0 + u 01 (y — yy) + u 02 (y — yy) 2 ]i,y is a piecewise continuous function from the wall to 
the free stream (See equation A. 18). 

Figure 6 compares u with the Blasius solution at the airfoil leading edge (x = .01). 
The discrete velocity clearly shows the effect of the strong leading edge singularity. The 
Blasius similarity solution, on the other hand, does not account for the singular behavior 
of the flow field in this region. It is apparent from these results that two cells are sufficient 
to resolve the leading edge boundary layer profile. 

In Figures 7 and 8 we compare with the Blasius solution at four different locations 
along the airfoil. The agreement steadily improves as x increases. This trend continues 
up to about x = .85, where trailing edge and wake effects become significant. This can be 
seen from the results shown in Figure 9 a., where we plot the maximum deviation of u 
from the Blasius solution as a function of x. 

In Figure 9 b., we show the improved accuracy that is obtained by refining the mesh 
spacing in the y direction. The exponential y stretching for the five grids used in the figure 
was 12.5%, 9.5%, 7.1%, 5.3%, and 3.8%, respectively. The maximum deviation of u from 
the Blasius solution is less than 5 x 10“ 3 over a significant portion of the airfoil for all of 
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the grids except the coarsest one. 

Jn Figures 6-8 we compared u with the Blasius solution at fixed nodal points x$. 
This enabled us to present a streamwise velocity profile which is continuous across cell 
interfaces. However, it only required using the three discrete Taylor coefficients u 0 0 , u 01 , 

and u 0 , 2 . In Figures 10 a. and b., we compare u(x,- + fj j 1 -) with the Blasius 

solution (i. e., we compare at the lower right hand corner of each solution element). This 
comparison requires the use of all the discrete Taylor coefficients associated with u. The 
mesh sizes associated with Figures 10 a. and b. are indicated below the figures. The 81 
x 22 grid is the same as the one used in Figure 9 b. (Ax = .02, A y w = .0015 with 9.5% 
exponential y stretching). The 110 x 28 grid is refined in both the x and y directions 
(Ax t .e. = -007, 2.5% exponential x stretching away from the trailing edge, A y w = .0015, 
with 3.8% exponential y stretching). The grids are intended to be “coarse” and “fine”, 
respectively. 

The accuracy demonstrated in Figures 10 a. and b. is comparable to that shown in 
Figures 7 and 8. This suggests that u uniformly approximates the exact solution u over 
the entire solution element. Further evidence for this is provided by the results shown in 
Figures 10 c. - f. In Figures 10 c. and d., comparison is made with the Blasius solution 
at the cell center (i. e., u 0 0 is compared with the Blasius solution). Figures 10 e. and f. 
show the “error” (i. e., the deviation from the Blasius solution) of the numerical results in 
Figures 10 a. - d. 

The above results confirm that u does indeed uniformly approximate the exact solution 
u throughout the solution element. This is a direct consequence of requiring the discrete 
flux vectors to identically satisfy the governing conservation law. Since the conservation 
law is satisfied identically, and not just at a point, there is no preferred location within a 
solution element. Thus, the accuracy of the discrete approximation is essentially uniform. 

One final set of results is presented in Figure 10 g. Here we compare the value of u 
from adjacent solution elements at a common point on their interface. The maximum 
difference in the value of u corresponding to the i = I and i = / + 1 locations is less than 
2.24 x 10 -3 . These results correspond to the “coarse” grid. For the “fine” grid, a similar 
calculation showed a maximum difference which is less than 7.7 x 10 -4 . Thus, u is nearly 
continuous across cell interfaces. 

We conclude our discussion with a presentation of numerical results from the trailing 
edge region. Our motivation was to determine how fine a mesh spacing is required to 
accurately resolve the trailing edge pressure singularity. The grids used for this purpose 
axe shown in Figure 11. The y spacing was the same for all five grids. Each grid has 
Ay w = .0015, 9.5% exponential y stretching, and free-stream mesh boundaries located at 
±y = .027. The mesh in Figure 11 a. has uniform x spacing, while the other four grids 
are refined in the trailing edge region with exponential stretching. The value of Ax at 
the trailing edge was .02, .01, .009, .008, and .007, respectively. This is denoted below 
each figure by “dxte = .02”, etc. The exponential stretching away from the trailing edge 
was the same in the wake as on the airfoil. The downstream boundary was located at 
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approximately x = 1.5 for all cases. 

In Figures 12 and 13 we present the numerical results. Figure 12 shows the pressure 
coefficient as a function of x. The results correspond in order to the five grids shown in 
Figure 11. It is remarkable how little mesh refinement is required to resolve the pressure 
singularity. Our results indicate a minimum value of —.014 for the pressure coefficient. 
This agrees well with the recent calculations of Srinivasan and Rubin. 9 

Figure 13 shows the streamwise velocity profile at four different x locations in the 
trailing edge region. The results correspond to the last solution element on the airfoil and 
the first three solution elements in the wake. 

All of the calculations associated with the present work were performed on the Cray 
YMP at the NASA Lewis Research Center. The CPU times ranged from 30.3 seconds 
for the coarsest grid (81 x 16) to 170 seconds for the finest grid (110 x 28). Most of the 
calculations took about 80 CPU seconds. Each case was started from uniform flow and 
converged in nine or ten Newton iterations to a maximum residual error which was less 
than 10 -1 °. 

The above CPU times can be considerably reduced by using a previous solution, rather 
than uniform flow, as the initial guess for Newton's method. Further reductions can be 
achieved by using equations A. 9 — A. 11 to eliminate additional variables from the discrete 
system of equations that must be solved. If three more variables are eliminated, the block 
sizes associated with the Jacobian matrix can be reduced from 10 Nj to INj. Since the 
operation count for Gaussian elimination is 0(n 3 ), where n is the block size, the CPU 
times could be reduced by a factor of ( y-) 3 « 2.9. Combined with using an improved 
initial guess, we estimate that the CPU times can be reduced by a factor of about five 
without making any changes to the matrix solution technique. 


Summary 

In this paper we have presented a new numerical scheme for the solution of the thin 
airfoil boundary layer problem. The results presented above have shown (i) the ability 
of the scheme to accurately resolve the thin airfoil flow field on grids which are much 
coarser than those used by conventional numerical methods, (ii) the uniform accuracy of 
the discrete solution throughout the solution element, and (iii) a nearly continuous discrete 
solution across cell interfaces. 

In Section II, using both an integral and differential approach, we rederived the discrete 
flux conservation equations presented in [3]. Here we presented a simpler method for 
deriving the flux expressions at cell interfaces, and provided a systematic and rigorous 
derivation of the conditions used to simulate the differential form of the conservation laws. 
A generalized concept of convergence was introduced, and necessaiy conditions for the 
order- JV convergence of a discrete approximation were derived. In addition, necessary and 
sufficient conditions for a discrete approximation to satisfy a conservation law in E 2 were 
presented. An ideal error bound on the discrete solution was also derived. 
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We conclude with the following remarks. First, the theoretical results established in 
Section II C. axe applicable to any conservation law in E 2 , and their analogues in higher 
dimensions readily follow. Second, extension of the present scheme to three dimensions 
is straightforward and follows naturally from the 2-D formulation. Third, higher-order 
schemes for the Navier-Stokes equations may be very efficiently constructed using the 
methodology and theory presented in Section II. Fourth, a faithful simulation of the con- 
servation laws associated with the Navier-Stokes equations requires the rigorous enforce- 
ment of both local and global flux conservation. This is accomplished most naturally 
through an integral formulation. At the same time, the integral and differential forms of 
the conservation laws cannot be divorced from each other. Thus, a fundamental tenet of 
the space-time solution element method is the rigorous enforcement of both the integral 
and differential forms of the governing conservation laws. Finally, the results presented 
in Section IV have clearly demonstrated the convergence of the discrete solution as the 
mesh spacing is refined. The uniform accuracy over the solution element suggests that the 
convergence is of order two. An analysis of the order of convergence will be provided in a 
forthcoming paper. 
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Balance of Mass and Momementum Fluxes Across Vertical Interfaces: 
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Balance of Mass and Momementum Fluxes Across Horizontal Interfaces: 
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Continuity of u Across Horizontal Interfaces: 
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Airfoil Boundary Conditions: 
Lower Surface 
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Upstream Boundary Conditions: 
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Downstream Boundary Condition: 
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Figure 3. Computational grid for a thin (flat-plate) airfoil. 
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Figure 5 a. 81 x 20 computational grid with 12.5% exponential y stretching. 


Figure 5 b. 81 x 28 computational grid with 3.8% exponential y stretching 
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Figure 6. Comparison of numerical results and Blasius solution at airfoil leading 
edge (x = .01). Calculations performed on an 81 x 20 grid with 12.5% exponential 
y stretching. 
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Figure 9 a. Maximum deviation of predicted streamwlse velocity from Blaslus solution 
as a function of x. 81 x 20 grid with 12.5% exponential y stretching. 



Figure 9 b. Effect of mesh refinement on deviation from Blaslus solution. 
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Figure 10 a. Comparison of numerical results and Blaslus solution at x = .82. 
Comparison Is at the lower right hand comer of each solution element. 81 x 22 grid. 
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Figure 10 b. Comparison of numerical results and Blaslus solution at x = .822. 
Comparison Is at the lower right hand comer of each solution element. 110x28 grid. 
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Figure 10 c. Comparison of numerical results and Blaslus solution at x = .81. 



Figure 10 cL Comparison of numerical results and Blasius solution at x = .81 6. 
Comparison is at the cell center of each solution element. 110 x 28 grid. 
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Figure 10 e. Difference between discrete u velocity and Blasius solution at cell 
center and lower right comer of solution elements, x = .81, .82. 81 x 22 grid. 



Figure 1 0 f. Difference between discrete u velocity and Blasius solution at ceil 
center and lower right comer of solution elements, x = .81 6, .822. 110x28 grid. 
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Figure 1 0 g. Comparison of numerical results from adjacent solution elements 
and the Blasius solution, x * .82, 81 x 22 grid. 


Figure 11 a. 81 x22 computational grid with uniform x spacing, dxte = .02. 


Figure 11 b. 98 x 22 computational grid with 2% exponential x stretching. 
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Figure 11 c. 102 x 22 computational grid with 2.1% exponential x stretching, 
dxte = .009 
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Figure 11 d. 106 x 22 computational grid with 2.3% exponential x stretching, 
dxte = .008 



Figure lie. 110x22 computational grid with 2.5% exponential x stretching, 
dxte = .007 
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Figure 12 a. Pressure coefficient in the trailing edge region. 81 x 22 grid. 



Figure 12 b. Pressure coefficient in the trailing edge region. 98 x 22 grid. 
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Figure 12 d. Pressure coefficient In the trailing edge region. 106 x 22 grid. 
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Figure 12 e. Pressure coefficient in the trailing edge region- 1 10 x 22 grid. 



Figure 13 a. Predicted stream wise velocity profile In the trailing edge region 
from an 81 x 22 grid with uniform x spacing, dxte = .02. 







Figure 13 b. Predicted streamwise velocity profile In the trailing edge region 
from a 110 x 22 grid with exponential x stretching, dxte = .007. 
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